Einführung

Was ist R

R ist eine öffentliche Statistik-Programmiersprache, welche sich für eine breite Anwendung einsetzen lässt.

Was ist RStudio

RStudio ist in R-Kreisen die bekannteste Programmier-Entwicklungsoberfläche. Sie wurde von der Firma RSTudio entwickelt und es gibt eine open-saurce Version.

Warum R

  • open source Software
  • sehr gut dokumentiert, da Verwendung weit verbreitet
  • wiederkehrende Datenaufbereitungs-Prozesse lassen sich leicht automatisieren
  • sehr stark in data wrangling
  • unterschiedlichste Tools, die es dem Datenanalysten vereinfachen ganze Prozessketten selbst zu erstellen ohne andere Programmiersprachen zu können

GIS mit R

Libraries die auf GDAL zugreiffen

Vektor Libraries:

  • sf
  • sp (durch sf abgelöst)

raster libraries:

libraries für Visualisierungen:

Buch zur Einführung:

Ausgangslage

Die katholischen Kirche möchte wissen wie gross der Anteil der ständigen Wohnbevökerung an über 80 Jährigen in den Kirchgemeinden ist.

Inputdaten: Kirchengemeinde (shp); Stand 2019 städnigen Wohnbervölkerung (sph); Stand 2020

Methodik

  1. Kirchengemeinde: Stand 2019 => Stand 2020:

    1. Die Kirchgemeinde “Hirzel-Schönenberg-Hütten” wurde getrennt:
      => Hirzel hat sich der Kirchgemeinde Horgen angeschlossen.

      => Schöneberg und Hütten haben sich der Kirchgemeinde Wädenswil angeschlossen.

  2. Die Daten der Bevölkerungsstatistik werden auf die neuen Kirchgemeinden aggregiert.

Hinweis:

Da die Bevölkerungsstatistik auf eine Hektare aggregiert werden, stimmen die aggregierten Zahlen nicht ganz mit der realität überein. Wir verwenden die Bevölkerungsstatistik da die Daten OGD sind und ihr so das Beispiel durchspielen könnt. Wir haben die Daten eigentlich auf EGID genau und können auf dieser Basis auswertungen machen.

Berechnung

Vorbereiten des Setups

Es werden alle packages geladen, welche für die Aufbereitung verwendet werden. Zudem wird eine Funktion erstellt, welche die Shapefiles auf nicht valide Geometrien untersucht.

# Laden aller benötigten Packages
library(sf)         # Vektor-Library
library(dplyr)      # Library zur Datenaufbereitung
library(stringr)    # Library für regex
library(ggplot2)    # Library zur Plot-Erstellung
library(ggspatial)  # Library mit Map-Features (Nordpfeil usw.)
library(mapview)    # Library für interaktive Karten

# Hilfsfunktion um die Validität der Geometrien zu prüfen. Falls eine Geometrie kaputt ist, gibt es eine Warnung. Die Topologie kann jedoch nicht geprüft werden.
check_validity <- function(sf_df){
  if(FALSE %in% sf::st_is_valid(sf_df)){
    warning("There is an invalid geometry in the Dataset. Use st_make_valid() to repair the geometry.")
  }else{
    print("The geometry is valid.")
  }
}

Importieren der Datensätze

Der Import der Datensätze, in unserem Fall shapefiles, kann mit st_read() gemacht werden. Es können aber auch zum Beispiel gpkg, geojson und Daten aus einer DB eingelesen werden. Daten können auch, mit hilfe eines SQL Queries, partiell eingelesen werden.

# Importieren des Kirchgemeinden-Shapefiles des Jahres 2019
kg <- sf::st_read(paste0("../shapefiles/kirchgemeinden_2019.shp"), stringsAsFactors = FALSE, crs=2056)
## Reading layer `kirchgemeinden_2019' from data source `/home/b105ptk@ji.ktzh.ch/git/GISwithR/shapefiles/kirchgemeinden_2019.shp' using driver `ESRI Shapefile'
## Simple feature collection with 62 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 2669245 ymin: 1223895 xmax: 2716901 ymax: 1283343
## Projected CRS: CH1903+ / LV95
# Importieren des Gemeindegrenzen-Shapefiles des Jahres 2020
gem <- sf::st_read(paste0("../shapefiles/gemeinden_2020.shp"), stringsAsFactors = FALSE, crs=2056, query = "SELECT GDE_ID, ART_N FROM \"gemeinden_2020\" ")
## Reading layer `gemeinden_2020' from data source `/home/b105ptk@ji.ktzh.ch/git/GISwithR/shapefiles/gemeinden_2020.shp' using driver `ESRI Shapefile'
## Simple feature collection with 162 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 2669245 ymin: 1223896 xmax: 2716900 ymax: 1283343
## Projected CRS: CH1903+ / LV95

Schauen wir uns die Daten an

Das sf package hat eine eigene plot-Funktion welche fähig ist Karten zu plotten. Diese ist sehr einfach, es hilft aber um sich einen schnellen Überblich zu verschaffen. Mit mapview() gibt es eine Funktion, die erlaubt, eine Karte interaktiv anzuschauen. Dazu wird im backend leaflet verwendet.

# statischer Plot
plot(kg$geometry)

# dynamischer Plot
mapview(kg)

Überprüfen ob die Geometrien der beiden Objekten ok sind

Das sf-package bietet Funktionen an um die Geometrie jedes polygons zu prüfen und zu korrigieren.

check_validity(kg)
## [1] "The geometry is valid."
check_validity(gem)
## [1] "The geometry is valid."

Erstellen der Kirchgemeindegrenzen für das Jahr 2020

Extrahieren der benötigten Polygone

Für das Splitting der Kirchgemeinde, werden die politischen Gemeindegrenzen verwendet. Also der Teil der Kirchgemeinde, welcher in der politischen Gemeinde Horgen liegt, wird der Kirchgemeinde Horgen zugewiesen, der andere Teil der Kirchgemeinde Wädenswil. Um den Prozess relativ simple zu halten, extrahiere ich nur die benötigten Polygone.

# Filtern der Kirchgemeinde mit der Nummer 19 (Hirzel-Schönenberg-Hütten)
kg_19 <- kg %>% filter(Kirchnumme == 19)

# Filtern der politischen Gemeinden (Horgen, Wädenswil)
gem_h_w <- gem %>% filter(GDE_ID %in% c(293, 295))
# interaktive Darstellung
mapview(gem_h_w, zcol = "ART_N") + mapview(kg_19, zcol = "ART_TEXT", col.regions = c("red"))

Erstellen einer Mapping-Tabelle

Um dann die gesplitteten Polygone den richtigen Kirchgemeinden zuordnen zu können, wird eine Mapping-Tabelle erstellt.

# Erstellung einer Tabelle mit tibble
kg_fus_tibble <- tibble(
  GDE_ID = c(293, 295),
  Kirchnumme = c(45, 21),
  GEMEINDENA = c("Wädenswil", "Horgen"),
  ART_TEXT = c("Kirchgemeinde", "Kirchgemeinde"),
  ART_CODE = c(1,1)
)

# Darstellen der Tabelle
knitr::kable(kg_fus_tibble)
GDE_ID Kirchnumme GEMEINDENA ART_TEXT ART_CODE
293 45 Wädenswil Kirchgemeinde 1
295 21 Horgen Kirchgemeinde 1

Verschneiden der Kirchgemeinde mit den politischen Gemeinden

Mit dem Verschnitt der beiden Layers wird die Kirchgemeinde in mehrere Polygone unterteilt, welche die ID der darunterliegenden politischen Gemeinden enthalten.

# Verschneiden der Kirchgemeinde mit den politischen Gemeinde um die Kirchgemeinde zu trennen. 
#Der Output wird als Multipolygon ausgegeben
union_gem <- st_intersection(kg_19, gem_h_w) %>% st_cast("MULTIPOLYGON")

mapview(union_gem, zcol = "GDE_ID")
# Multipolygon -> single Polygon
union_single <- st_cast(union_gem, "POLYGON")

mapview(union_single, zcol = "GDE_ID")

Zusammenfügen der kleinen Sliver-Polygone mit den grossen Polygonen

Es gibt einige Sliverpolygone, welche den zwei grossen Polygonen zugeordnet werden.

# Die sliverpolygon kleiner als 1000m2 werden der Kirchgemeinde Horgen zugeordnet.
# Wir können das so machen, da es nur sliver-polygone an der Grenze zwischen der 
# Kirchgemeinde Horgen und der politischen Gemeinde Wädenswil auftreten.
union_merged <- union_single %>% mutate(area = as.numeric(st_area(.))) %>% 
  mutate(GDE_ID = ifelse(area < 1000, 295, GDE_ID)) %>% 
  group_by(GDE_ID) %>% 
  summarize()

mapview(union_merged, zcol = "GDE_ID")

Zuordnen der neuen Kirchgemeinden-ID’s zu den neu entstandenen Polygonen

Die neuen Polygone enthalten zur Zeit nur die politische Gemeinde-ID und die alte Kirchennummer. Mit der Mapping-Tabelle werden nun die neuen Kirchgemeinde-ID hinzugefügt.

# Die mapping Tabelle wird nun an die beiden neuen Polygone hinzugefügt um die Kirchennummer den Polygone zuzuordnen.
new_kirchid <- union_merged %>% 
  left_join(kg_fus_tibble, by = "GDE_ID") %>% 
  select(-GDE_ID)

Ersetzten der alten Kirchgemeinden mit den neuen Polygonen und mergen der Kirchgemeinden

Wir ersetzten nun die Kirchgemeinde mit der Nummer 19 durch die zwei neuen Polygone. Die neuen Polygone werden dann mit den bestehenden Kirchgemeinden gemerged um einheitliche Polygone zu erhalten.

# Ersetzen der alten Kirchgemeinde durch die zwei neuen Polygone.
# Dann werden die Kirchgemeinden zusammengefügt
kg_2020 <- kg %>% 
  filter(Kirchnumme != 19 | is.na(Kirchnumme)) %>% 
  bind_rows(new_kirchid) %>% 
  st_cast("MULTIPOLYGON") %>% 
  mutate(STICHTAG = "2021-01-01") %>% 
  group_by(ART_TEXT, ART_CODE, STICHTAG, GEMEINDENA, Kirchnumme) %>% 
  summarize()

check_validity(kg_2020)
## [1] "The geometry is valid."
plot(kg$geometry)
plot(kg_2020$geometry)

So, jetzt haben wir den Kirchgemeindestand vom Jahr 2020.

Laden der Punktdaten der Bevölkerungsstatistik

Um nun den Anteil der über 80 Jährigen zu berechnen, verwenden wir den OGD-Datensatz der Bevölkerungsstatistik. Da dies Hektardaten sind, ist es lediglich eine Annäherung. In unserer DB hätten wir Gbäude-Spezifische Daten, um das Beispiel nachspielen zu können, wird hier aber der OGD-Datensatz verwendet.

Da bereits die Anteile auf Hektarraster-Ebene berechnet wurden, müssen wir zuerst die Prozente wieder in Absolute Zahlen umrechnen.

# laden der Bevölkerungsstatistik
bevstatistik <- sf::st_read("../shapefiles/bevstatistik_punkte.shp")
## Reading layer `bevstatistik_punkte' from data source `/home/b105ptk@ji.ktzh.ch/git/GISwithR/shapefiles/bevstatistik_punkte.shp' using driver `ESRI Shapefile'
## Simple feature collection with 36256 features and 23 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 2669650 ymin: 1224850 xmax: 2716150 ymax: 1283350
## Projected CRS: CH1903+ / LV95
# Vorbereiten der Daten
bevstatistik_abs <- bevstatistik %>% 
  mutate_at(vars(matches("_P")), ~round(PERS_N/100*.)) %>% 
  rename_at(vars(matches("_P")), ~str_replace_all(., "_P", "_A")) %>% 
  mutate_at(vars(-geometry), ~ifelse(. == 9980 | . == -999, NA, .))

Zuordnen der Punktdaten zu den Kirchgemeinden

Die Punktdaten und Kirchgemeinden werden verschnitten um dann über die Kirchgemeinde aggregieren zu können.

# Verschneiden der Punktdaten mit den Kirchgemeinden
bev_within_kg <- sf::st_intersection(bevstatistik_abs, kg_2020)

mapview(bev_within_kg %>% filter(Kirchnumme==45))

Aggregieren der Daten und berechnen des Anteils der über 80 Jährigen an der Bevölkerung in den Kirchgemeinden

Die Punktdaten besitzen nun die Information der Kirchgemeinde und wir können die Daten nach dieser Variable aggregieren. Dann müssen wir noch den Anteil der über 80 Jährigen berechnen.

# aggregieren der Daten
anteil_ue80_pro_kg <- bev_within_kg %>% 
  st_drop_geometry() %>% 
  filter(!is.na(Kirchnumme)) %>% 
  group_by(Kirchnumme) %>% 
  summarize(anzahl_ue_80 = sum(J_80PLUS_A, na.rm = T),
            anzahl_tot = sum(PERS_N, na.rm = T)) %>% 
  ungroup() %>% 
  mutate(anteil_ue_80 = round(anzahl_ue_80/anzahl_tot*100, 1))


knitr::kable(head(anteil_ue80_pro_kg, 10))
Kirchnumme anzahl_ue_80 anzahl_tot anteil_ue_80
1 1136 19016 6.0
2 1143 26075 4.4
3 1148 22712 5.1
4 524 11770 4.5
5 795 13242 6.0
6 600 14330 4.2
7 1586 35178 4.5
8 1539 42131 3.7
9 1370 28268 4.8
10 2236 43619 5.1

Joinen der Aggregierten Daten an das Kirchgemeinde-Schapefile und erstellen einer Karte

Die Aggregierten Daten können nun dem Kirchgemeindenshapefile angefügt werden und wir können nun eine Choroplethenkarte erstellen.

kg_2020_with_anteil <- kg_2020 %>% 
  left_join(anteil_ue80_pro_kg, by = "Kirchnumme") %>% 
  select(Kirchnumme, anteil_ue_80)



# plot der Daten
ggplot(kg_2020_with_anteil) +
  geom_sf(aes(fill = anteil_ue_80), color = "white") +
  theme_bw()  +
  theme(panel.grid.major = element_line(colour = "transparent")) +
  scale_fill_continuous(name = "Anteil ü80 [%]") +
  theme(legend.position = c(0.87, 0.85)) +
  theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
        legend.key.height = unit(0.5, 'cm'), #change legend key height
        legend.key.width = unit(0.5, 'cm'), #change legend key width
        legend.title = element_text(size=9), #change legend title font size
        legend.text = element_text(size=7)) +
  coord_sf(datum = NA) +
  annotation_scale(location = "bl", width_hint = 0.5) +
  annotation_north_arrow(location = "bl", which_north = "true", 
        height = unit(1, "cm"), width = unit(1, "cm"),
        pad_x = unit(0.1, "in"), pad_y = unit(4.2, "in"),
        style = north_arrow_fancy_orienteering)